Pore Stabilization in Cohesive Granular Systems 
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Cohesive powders tend to form porous aggregates which can be compacted by applying an external 
pressure. This process is modelled using the Contact Dynamics method supplemented with a cohe- 
sion law and rolling friction. Starting with ballistic deposits of varying density, we investigate how 
the porosity of the compacted sample depends on the cohesion strength and the friction coefficients. 
This allows to explain different pore stabilization mechanisms. The final porosity depends on the 
cohesion force scaled by the external pressure and on the lateral distance between branches of the 
ballistic deposit r cap t. Even if cohesion is switched off, pores can be stabilized by Coulomb friction 
alone. This effect is weak for round particles, as long as the friction coefficient is smaller than f . 
However, for nonspherical particles the effect is much stronger. 



PACS numbers: 02.70.Ns, 45.70.Cc, 62.25.+g 
I. INTRODUCTION 



Cohesion plays an important role in modern powder 
and nano technology. It determines the compaction and 
sintering behaviour as well as the mechanical properties 
like yield under shear stress Q An important tool to un- 
derstand the behaviour microscopically is the modelling 
on the particle scale, i.e. discrete element methods ||. 

So far most computational studies have neglegted co- 
hesion between particles. This is justified in dry systems 
on scales where the cohesive force is weak compared to 
the gravitational force on the particle, i.e. for sand and 
coarser material, which collapses under its own weight 
into a random dense packing. Powders already show co- 
hesion effects: With decreasing grain diameter cohesive 
forces lead to highly porous media. For particle diame- 
ters in the nanometer range the cohesive force becomes 
the dominant force, so that particles stick together upon 
first contact. 

When sintering nano-powders care has to be taken that 
the grains do not grow. This is achieved in so-called sin- 
ter forging under high pressure at relativly low temper- 
atures In this process the highly porous powder gets 
compacted. In order to simulate this compaction process 
we use the contact dynamics method, which in principle 
[4| allows perfectly rigid particles with Coulomb friction 
[5,|J without regularization of these force laws. Whereas 
in soft particle molecular dynamics there already exist co- 
hesion models [|7|-|l0|] , in contact dynamics cohesive mod- 
els are just at the beginning [|ll] [TJ). 

Here we present contact dynamics results for two di- 
mensional systems of round particles, which can be re- 
alized in experiments by aggregates of parallel cylinders 
|l4|| . Cohesion as well as rolling friction were included, 
which turns out to be crucial to describe the high poros- 
ity of nano powders. In three dimensions, torsion friction 
is needed in addition |lq] . Here we explain how the differ- 



ent contact properties contribute to stabilizing the pores. 



II. CONTACT DYNAMICS SIMULATON 

A. Without Cohesion and Rolling Friction 

It is believed that the physics of dry dense granular 
matter with many lasting contacts is determined by vol- 
ume exclusion and Coulomb friction. Both contact laws 
are nonsmooth, i.e. they do not determine the forces 
as functions of state variables (positions and velocities). 
The forces at a nonsliding contact are reaction forces in 
the sense that they have to compensate all external forces 
which would violate the constraints of volume exclusion 
and zero relative velocity. The contact dynamics method 
Plppl] is designed such that it determines these con- 
straint forces in every time step. Elastic deformations of 
the particles need not be taken into account, hence the 
choice of the simulation time step is not coupled to the 
stiffness of the particles as it is in soft particle molecular 
dynamics. 

For two particles in contact, the constraint force can be 
determined analytically, if the externally applied forces 
are known. But typical dense systems consist of many 
particles involving many contacts forming contact net- 
works, and calculating all contact forces is a global prob- 
lem, because the force at a contact influences the neigh- 
boring contacts and so on. To find a solution (which is 
not necessarily unique, cf. e.g. frlTf ), usually an iterative 
procedure is applied: The force at each contact is cal- 
culated in random order by considering the preliminary 
forces of adjacent contacts as already correct (which al- 
lows the same treatment as for external forces). This 
iteration goes on until a convergence criterion is fulfilled. 
A sensible choice is to demand the relative force changes 
at every contact to be below a given treshold during a 
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FIG. 1. Signorini graph (normal force Fn vs. distance d) 
for a cohesive contact: Distances above zero are allowed while 
interpenetration is avoided by means of _Fn > 0. In contrast to 
the usual Signorini graph, a negative force (attractive force) 
is allowed within the cohesion range dc- (Fat lines show the 
allowed values.) 
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FIG. 2. Coulomb graph (tangential force Ft vs. tangential 
velocity vt) for a cohesive contact: The tangential force Ft 
hinders the contact from sliding up to the Coulomb threshold, 
above the contact starts sliding. This threshold is increased 
due to the "normal" cohesion. (Fat lines show the allowed 
values.) 

certain number of consecutive iterations. 



B. Cohesion Model 

In this section the cohesion model we use ("normal 
cohesion" ) is described. For simplicity it has only two 
parameters, a constant attractive cohesion force Fc in 
normal direction, which acts over a short range dc de- 
termining the cohesion energy Ec = Fcdc (Fig- |l|). The 
particles are still considered as perfectly rigid: The co- 
hesive interaction does not deform them. A contact can 
only open, if an external pulling force exceeds the thresh- 
old Fc and performs work E larger than Ec, so that the 
particles separate with a kinetic energy E — Ec ■ While 
Radjai et al. [jllj use the same cohesion model, Jean et 
al. |12] proposed a different implementation, the FCR 



model, which has more parameters. 

The introduction of a force scale Fc into Signorini's 
condition (Fig. [I]) leads to a numerical complication. 
This is related to the existence of "shocks": Because of 
the perfect rigidity of the particles, a finite momentum 
Ap can be transmitted instantaneously, if the connected 
cluster of particles, to which the contact belongs, collides 
with some other particle or cluster. This corresponds to 
an infinite normal contact force in form of a Dirac-pulse 
Fjv = ApS(t). Due to the time discretization by steps 
of At, though, one gets a finite Fn = Ap/ At during the 
time-step containing the shock, which cannot be distin- 
guished from a force evolving continuously at a lasting, 
non-shocked contact. Whereas for Fc — shocks and 
persisting contacts can be treated in the same manner [^) , 
it depends on At, whether a contact with given Fc > 
and Ec > opens or not, if it is shocked (Ap < 0). 
For large At it may happen, that \Ap/At\ drops below 
the threshold Fc, so that the contact cannot open. In- 
deed we find that the number of contacts which open in 
a simulation run decreases with increasing At. 

However, this is not the only source of systematic er- 
rors related to the discrete time step: Probably most im- 
portant is that the iterative determination of the forces 
is only accurate within a given tolerance, so that the 
constraints are not perfectly obeyed. Hence a time step 
in general leads to a small overlap between the parti- 
cles, which depends on At. The third systematic error 
occurs in time steps during which a contact opens. As 
the cohesion force is acting over the whole time step, al- 
though d = dc is reached some time in between, the 
cohesion energy dissipated during opening the contact is 
slightly overestimated. All three systematic errors be- 
come smaller for smaller time steps. 

As an example the final piston position as a result of 
a simulation of uniaxial consolidation by a fixed external 
force is shown in Fig.|| as a function of the time step. For 
decreasing time step the simulation result systematically 
increases. In our simulations in the next sections we use 
a time step At = 0.01, which leads to an average overlap 
between the particles of about 5% of their radius. This 
explains why in Fig.|^ the final piston position for At — 
0.01 is about 5% smaller than the extrapolated value at 
Deltat = 0. To get a lower percentage one would have 
to use smaller time steps which leads to an increase of 
computation time. We only compare simulation results 
which were obtained with the same time step. 

We also did quasistatic simulations, where shocks are 
eliminated by resetting the velocities to zero after every 
time step. Then the results are different from the case 
of finite compaction speed: One gets more porous struc- 
tures within the compaction process by a piston with 
an external load. Due to neglecting inertia effects the 
system's reaction force onto the piston is below the ex- 
ternal force on the piston at any time. In the "end" both 
forces are almost equal, but a low difference is left lead- 
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FIG. 3. The final position of the compacting piston after 
consolidation is influenced by the simulation time step used. 



ing to a sustained movement of the piston. Thus, the 
final static configuration as found in the simulations pre- 
sented in this paper, cannot be reached via a quasistatic 
simulation in finite computation time. This simulation 
methods would be applicable for simulating systems un- 
der constant strain only. 

The opening of a contact needs usually several time 
steps, in which the pulling force exceeds Fq- The pulling 
force minus Fq has to perform the work Eq necessary 
to reach the distance d — dc, where the contact breaks. 
In our model a contact which started to open, but at 
which d — dc has not yet been reached, is not pulled 
back by the cohesive force, if the pulling force becomes 
weaker than Fq (concerning nano-particles a sinter neck 
which was pulled by an external force to a thinner neck is 
not built up to its former width again without sintering) . 
Such a weakened but not yet broken contact can only be 
strengthened again (decrease of d), if the particles are 
pushed together externally. 

For comparison we also did some simulations with a co- 
hesion model, in which the cohesion force is able to pull 
the particles together again, if the contact had not been 
broken. A contact which has been slightly opened by an 
external force (i.e. forces from the other contacts on the 
particles) is pulled back to zero gap between the parti- 
cles again (the normal force not exceeding the threshold 
—Fq). As a result a contact, which did not break, has 
a chance to relax, instead of remaining weakened. Sim- 
ulating this model gives similar results as presented in 
this paper. The main difference is due to the roughness 
of the upper envelope of the ballistic deposits which we 
choose as initial configuration for the simulation of the 
compaction process (Fig. ||). For high cohesion values 
the system really stays nearly unchanged, whereas in the 
simulations presented in the following the system com- 
pacts, until the upper envelope has become flat. 

Cohesion does not only counteract the opening of a 
contact but also hinders it to become sliding. Therefore 
the Coulomb friction law needs to be modified. In order 
to avoid introducing additional parameters, we assume 
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FIG. 4. Rolling condition: Similar as in the Coulomb graph 
the local torque Ti oc hinders the two particles from rolling 
against each other up to a threshold. Above this threshold 
the contact becomes rolling. 



that the Coulomb friction is simply raised to h(F_n +Fc) 
like in |ll| (Fig. |J). An alternative modification would 
be to increase the friction coefficient /i depending on Fe- 
in Sec. Ill we show, that this by itself would lead to the 
stabilization of pores even if "normal cohesion" would be 
absent. 



C. Rolling Friction Model 



In the case of round particles, rolling of two particles 
on each other must be considered. Its special importance 
in the context of cohesive contacts (as opposed to purley 
repulsive particles) can be explained as follows: Cohe- 
sion leads to stabilization of chains (i.e. particles with a 
coordination number of two), and thus allows for a less 
compact structure. The latter effect, however, is possi- 
ble only if rolling is suppressed, otherwise the chains are 
floppy and will fold. 

Introducing rolling friction to the model means to al- 
low for a local torque. But in principle, the point contact 
formed by two adjacent rigid discs (or spheres) cannot 
exert a torque. This may be compared to the fact that 
the microscopic origin of Coulomb's law of friction is not 
obvious for a contact with zero area. In both cases we re- 
gard the dimensions of the contact area and deformation 
zone as negligibly small compared to the particles. 

The condition for rolling is similar to Coulomb's law 
for sliding (Fig. ||): The contact can bear a local torque 
up to the threshold (i T (F N + F c ) (cf. Fig. g). If this is 
exceeded, the contact becomes a rolling one and hence 
exhibts a local relative angular velocity w T . According 
to experimental results |Q8[, the coefficient for rolling 
friction ji x is not chosen to be velocity dependent as it 
would be for a viscoelastic material |lj] , but is regarded 
as a phcnomcnological material constant which includes 
a (microscopical) length scale. 
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FIG. 5. Initial configuration: identical sperical particles are 
ballistically deposited. 

For the simulations shown in the following, we con- 
sider only discs of identical size and mass. Only then the 
equations of motion for the relative angular velocity w r 
on the one hand, and for the relative translational ve- 
locity on the other, are not coupled, which simplifies the 
implementation of rolling friction significantly. 



III. SIMULATIONS RESULTS 

A. Initial Configuration 

In principle there are two general methods to pro- 
duce nano-particles: the first is to fracture bulk mate- 
rial. The second one builds nano-particles by chemical 
precipitation or condensation in a liquid or gas environ- 
ment. While the production in large quantities is usually 
done in liquids, the production in the gas phase has some 
advantages (higher pureness of the material and sharper 
grain size distribution). To extract nanoparticles out of 
a gas flow one can use a filter where the particles are 
deposited on. Simulations of the filter process by a fiber 
network show a finger-like structure of the deposit |^0|,|lJ . 
Similar structures are obtained on a flat surface by ballis- 
tic deposition: Particles fall from the top randomly and 
stick to the first particle they reach within a capture ra- 
dius r capt around the particle. Such a ballistic deposit is 
shown in figure || for two dimensions j2^] . These struc- 
tures are not fractal, i.e. the porosity is not depending 
on the system size but on the capture radius. The cap- 
ture radius is a measure for the average distance between 
the branches so that the volume fraction of a ballistic de- 
posit increases about linearly with the reciprocal of the 
capture radius (Fig. ^) . For monodisperse spherical par- 
ticles the minimum capture radius is two particle radii, 
i.e. the minimal distance between the centers of mass two 
particles can have. In the following sections different bal- 
listic deposits with periodic boundary conditions in lat- 
eral direction are compacted by a piston with constant 
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FIG. 6. The density (volume fraction) of the ballistic de- 
posited configurations decreases with increasing capture ra- 
dius (here plotted against l/r capt ). 

load in vertical direction. 



B. Influence of different contact properties 

Initial configuration is a ballistic deposit with capture 
radius of 2.5 particle radii (Fig. ||) containing 433 parti- 
cles. After the compaction process by applying an exter- 
nal force on the upper piston one finally reaches a final 
equilibrium state. Figure |?j shows the final state for differ- 
ent contact properties used within the simulation. If one 
uses Coulomb friction (C) only (in addition to excluded 
volume interaction) one ends up with a high volume frac- 
tion (78%). Using rolling friction in addition (CR) a vol- 
ume fraction of 77% is reached. Compared to a perfect 
triangular lattice with volume fraction 7r/\/l2 « 90.7% 
there are only a few larger pores in these two configura- 
tions. This can be explained by the fact that force chains 
must be stabilized by adjacent particles avoiding the lat- 
eral movement of the chains |25| ]. Coulomb friction and 
cohesion (CC) lead to larger pores (see Fig. 0) and thus 
a volume fraction of 68%. The stabilization mechanism 
of force chains (and thus pore stabilization) is a different 
one: Force chains of compressive forces (gray lines) are 
stabilized by nearby force chains of tensile forces (black 
lines) so that the piston load is carried by the system. 
The higher thickness of the force lines indicate that the 
stresses in the system are much higher than in the ab- 
sence of cohesion. However, unhindered rolling leads to a 
destabilization of single particle chains. Therefore addi- 
tional rolling friction (CRC) leads to the highest poros- 
ity (volume fraction of 51%) because the system includes 
stable single particle chains. In this case there is sta- 
bilization of three degrees of freedom: the separation of 
the particles, the lateral movement as well as the rolling 
against each other. These three degrees of freedom are 
stabilized by the combination of cohesion, Coulomb fric- 



4 




FIG. 7. Configurations after compaction with different con- 
tact properties: Coulomb friction only (C); Coulomb and 
rolling friction (CR); Coulomb friction and cohesion (CC); 
Coulomb and rolling friction and cohesion (CRC); gray and 
black lines show compressive and tensile forces, respectively. 
The line thickness is a measure for the force value. The size 
of the black dots at the contacts indicates the magnitude of 
the opposite torques exerted by the contact on the adjacent 
particles. 



FIG. 8. Ballistic deposit before and after compaction by 
a constant external force on the piston. Contact properties 
include Coulomb friction, rolling friction and cohesion. 



tion, and rolling friction. 



C. Compaction process 



In order to get insight into the compaction process 
itself the temporal behaviour of the system is studied 
in this section. We simulated a ballistic deposit with 
capture radius r capt = 2.5 containing 2746 monodis- 
perse spherical particles (Fig. ||) with periodic bound- 
ary conditions in lateral direction. We compared the 
compaction dynamics for two different contact proper- 
ties between the particles: coulomb friction (C) only 
leads to the most compact final configuration, whereas 
Columb, rolling friction and cohesion (C RC) l ead to the 
most porous final configuration (see sec. Ill B ) . In both 
cases the time evolution of the piston position (Fig. ||) 
shows three phases of the compaction process. In the 
first phase the piston accelerates. Initially it is in contact 
with only one branch of the ballistic deposit. Only later it 
sweeps up substantial mass. Then, in the second phase 
the piston moves with a nearly constant velocity. The 
momentum transfered to the system by the force acting 
on the piston (in the cohesive case partly compensated 
by a reaction force of the opposite side of the container) 
is consumed by a linearly increasing mass swept up by 
the piston: 



F, 



piston 



reaction V 



(i) 



In the third phase the piston finally reaches a constant 
position (only neglectable small oscilations due to quasi- 
elasticity [j|): F pis ton = ^reaction- The final position of 
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FIG. 9. Time evolution of the piston's position during the 
compaction process for Coulomb friction only (C, dashed) and 
Coulomb and rolling friction and cohesion (CRC, full line) in 
comparison. 
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FIG. 11. Scaling plot for different cohesion strength shows 
alignment for Systems with different capture radius for small 
cohesion pressure ratios only (cp. Fig. |l0| ). For strong cohe- 
sion the final densities are mainly given by the initial densities. 



i10"' r 




1(T r 



force Fc by the pressure exerted by the piston, F p i S t on /L x 
in the two dimensional case. In the beginning an (almost 
linear) increase can be seen. For large cohesive force the 
scaled final piston position reaches a constant value. In 
this region there is almost no compaction. The fact that 
FcL x / Fpiston is the relevant quantity for the compaction 
process implies that the typical distance of strong force 
lines is not depending on the system size, as is also found 
for systems without cohesion |23j|. 



E. Scaling behaviour for different initial densities 



FIG. 10. Scaling plot for different cohesion strength shows 
alignment for different Systems. The systems are different in 
size and aspect ratio but have the same capture radius for 
ballistic deposition and thus approximately the same initial 
density. 



the piston is differ ent du e to the different stabilization 
mechanism (see sec. [II B ). 



D. Influence of cohesion strength 



The final piston position of different systems compact- 
ified by a constant external force is analyzed. At the 
contacts Coulomb and rolling friction, as well as cohesion 
with different values of Fc are acting. The initial systems 
are ballistic deposits with capture radius r capt = 2.5 and 
periodic boundary conditions. They consist of different 
numbers of particles and have different aspect ratios. We 
found that all the data collapse onto a single curve, Fig- 
ure M, if the final piston position is scaled by the final 
piston position without cohesion, y m i n , and the cohesive 



In the previous section the compaction of ballistic de- 
posits with the same capture radius was investigated. 
Now three ballistic deposits of the same size in vertical 
and lateral directions, but different capture radius are 
compacted by a piston with a fixed external load. Thus 
the initial systems contain different numbers of spherical 
particles, namely 2777 (r capt = 2), 2111 (r capt = 3) and 
1624 particles (r cap t = 4). The initial densities scale with 
lAcapt (see fig. ||). The same plot as in figure [To| (each 
curve averaged over 10 runs) leads to a data collapse for 
small ratios of cohesion and external pressure on the pis- 
ton only. In this region the different initial densities do 
not play an important role. For high cohesion values the 
final densities remain essentially the initial densities, so 
that there is no scaling in this region. 



F. Pore stabilization without cohesion 

The results presented in Figs. [To|and pT| suggest the fol- 
lowing physical picture: The external load of the piston 
must be carried by a set of strong force lines, which have 
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FIG. 12. Initial and final configuration after compaction 
with friction coefficient /j, = 3.8 and rolling friction coefficient 
fi r = 2.0 instead of cohesion. 
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FIG. 13. Increased friction coefficient p leads to pore sta- 
bilization and thus to lower density. Additional rolling friction 
amplifies this effect. 



a typical distance. The external load per force line is pro- 
portional to the pressure F p i ston /L x . It can be viewed as 
the destabilizing force along a force line and must be bal- 
anced by the stabilizing influence of another force, which 
in the previous sections was related to cohesion (and the 
fixed rolling friction), so that the scaling variable was 
ix-fc/^piston- This variable had to be big enough to pre- 
vent compaction. 

The question we want to address in this section is, 
whether friction forces alone can provide the stabiliza- 
tion as well. It is plausible to identify the external load 
per force line with the typical normal force at a contact 
along the force line, F n ~ F p i a t on / L x . The ratio between 
the stabilizing friction force Ft = fiF a and the destabiliz- 
ing force would then be L K F t /.Fpistou ~ M- This argument 
suggests that strong enough Coulomb friction may stabi- 
lize pores. This is indeed the case, as Fig. fi"3| shows. This 
effect is not very high (maximally about 10%), because 
the Coulomb friction cannot provide stabilization against 
buckling of the force lines. Therefore strong force lines 
need weak forces from the side as pointed out in [Q . Al- 
ternatively, rolling friction may stabilize the force lines 
against buckling. In combination with Coulomb friction 
this allows much larger pores, as the dashed curve in Fig. 
|l3| shows. However, the pore geometry is totally differ- 
ent from the one obtained in cohesive materials (compare 
Figs. | and |T^) . In the absence of cohesion large pores are 
found underneath arches. Presumably strong force lines 
stop the motion of grains on the upper side, while below 
pores open up due to inertia effects. Cohesion would lead 
to correlated motion of clusters of grains and would also 
prevent the inertial rupture of the structure underneath 
an arch. 

Without the use of rolling friction higher porosity will 
be reached by the use of non-sperical particles. To show 
this effect we simulated a a system of about 1700 parti- 
cles consisting of two different convex polygonal particle 
types and one spherical particle type. The diameter of 
the particles is chosen within the same range. The ini- 
tial configuration is a random loose packing, where no 
particles are in contact. After precompaction to a denser 
state the system is compacted by a piston with constant 
load using Coulomb friction coefficient — 0.5 (Fig. fuj ), 
respectively without Coulomb friction (Fig. |TJ). The two 
final configurations (Fig. [IJ, |l5|) have different porosity. 
It is higher, if Coulomb friction is present. Due to the 
shape disorder this effect is stronger than for the sys- 
tem of spherical particles. Here one ends up with a value 
y/l/min — 1 ~ 0.17 compared to the value y/y m in — 1 ~ 0.09 
at (i = 0.5 for the ballistic deposits of monodisperse sper- 
ical particles (Fig. |lj). One concludes that the corners 
of the particle have a similar effect as rolling friction for 
spherical particles. 

Thus, for cohesive non-sperical grains one expects a 
higher porosity depending on the shape of the particles, 
due to a similar stabilization mechanism as we found for 
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sperical particles using rolling friction. 



IV. SUMMARY. 



Different contact properties lead to different configura- 
tions after compaction. Using only Coulomb and rolling 
friction the final packing is an almost compact structure. 
The effect of cohesion on the porosity is higher if one 
uses Coulomb and rolling friction in addition. In that 
case single particle chains are stable and thus the pack- 
ing is highly porous. Important is the strength of the 
cohesion: Low cohesive forces lead to similar packings as 
without cohesion. Important is the relation between the 
stabilizing (cohesion) and destabilizing (pushing force) 
forces. Here one has to take into account that forces are 
extremly inhomogeneously distributed in the system so 
that one must devide the external pushing force by the 
average distance between strong force lines. This picture 
is confirmed by applying it to a system with increased 
values for the friction coefficient (above 1), which also 
leads to higher porosity. 
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